# This files estimates the discontinuity for variables expected to be unaffected to children's eligibilty
load("agg_day_year_all.rdata")
load("n_day_year_all.rdata")

## commented out for replication files 

#  install.packages("ggplot2")
#  install.packages("rmeta")
#  install.packages("reshape2")
#  install.packages("dplyr")
#  install.packages("rdrobust")

library(dplyr)
library(rmeta)
library(rdrobust)
library(ggplot2)
library(reshape2)

# write function for two samples proportion test

pr_mean_se <- function(var, by){
  p_1 <- var[by == 1]
  p_0 <- var[by == 0]
  mean_1 <- summary(p_1)[4] 
  mean_0 <- summary(p_0)[4]
  pr_mean <- mean_1 - mean_0
  pr_se   <- sqrt( mean_1*(1-mean_1) / length(p_1) + 
                     mean_0*(1-mean_0) / length(p_0) )
  p       <- 2*(1-pnorm(abs(pr_mean)/pr_se)) 
  n_under <- length(p_0)
  n_over  <- length(p_1)
  return(c(pr_mean, p, pr_se, mean_0, n_under, n_over))
}

# choose bandwidth based on when p < 0.05 parents with children 14 through 17 years or 19 through 22 years
ps <- data_frame(par.nat = rep(NA, 183),
                 par.fem = rep(NA, 183),
                 par.edu = rep(NA, 183),
                 chi.fem = rep(NA, 183),
                 chi.old = rep(NA, 183),
                 chi.coh = rep(NA, 183))



# estimate the effect at all bandwidths 
close_ests <- data_frame(par.nat = rep(NA, 183),
                         par.fem = rep(NA, 183),
                         par.edu = rep(NA, 183),
                         chi.fem = rep(NA, 183),
                         chi.old = rep(NA, 183),
                         chi.coh = rep(NA, 183))

close_se <- data_frame(par.nat = rep(NA, 183),
                       par.fem = rep(NA, 183),
                       par.edu = rep(NA, 183),
                       chi.fem = rep(NA, 183),
                       chi.old = rep(NA, 183),
                       chi.coh = rep(NA, 183))

close_tu <- data_frame(par.nat = rep(NA, 183),
                       par.fem = rep(NA, 183),
                       par.edu = rep(NA, 183),
                       chi.fem = rep(NA, 183),
                       chi.old = rep(NA, 183),
                       chi.coh = rep(NA, 183))

n_under  <- data_frame(par.nat = rep(NA, 183),
                       par.fem = rep(NA, 183),
                       par.edu = rep(NA, 183),
                       chi.fem = rep(NA, 183),
                       chi.old = rep(NA, 183),
                       chi.coh = rep(NA, 183))

n_over   <- data_frame(par.nat = rep(NA, 183),
                       par.fem = rep(NA, 183),
                       par.edu = rep(NA, 183),
                       chi.fem = rep(NA, 183),
                       chi.old = rep(NA, 183),
                       chi.coh = rep(NA, 183))

## Label mock variables

mockvars_labels <- c("Parents is native",
                     "Parent is female",
                     "Parent has bachelor's degree",
                     "Child is female",
                     "Child i oldest",
                     "Parent and child cohabit")

mock_vars <- c("par_native",
               "par_female",
               "par_degree",
               "chi_female",
               "chi_oldest",
               "par_cohabit")

data_master <- 
  left_join(data_day_year, data_n_year, by = c("days", "year")) 

for (j in 1:6){
  ## Create large dataset 
  var <- mock_vars[j]  
  data <- 
    data_master %>%
    ungroup() %>%
    mutate(par_vote  = unlist(data_master[,var])) %>% 
    select(year, days, par_vote, n)
  for(i in 1:183){
    data_close <- 
      data %>%
      filter(abs(days) < i) %>%
      mutate(voted     = round(par_vote*n),
             abstained = round((1-par_vote)*n),
             treated   = days > 0) %>% 
      select(days, treated, voted, abstained)
    
    data_voted <- 
      as.data.frame(lapply(data_close, 
                           function(x,p) rep(x,p), 
                           data_close[["voted"]])) %>%
      mutate(voted = 1)
    
    data_abstained <- 
      as.data.frame(lapply(data_close, 
                           function(x,p) rep(x,p), 
                           data_close[["abstained"]])) %>%
      mutate(voted = 0)
    
    data_close <- 
      rbind(data_voted, data_abstained)
    
    test <- with(data_close, pr_mean_se(voted, treated))
    
    close_ests[i,j] <- test[1]
    ps[i,j] <- test[2]
    close_se[i, j]  <- test[3]
    close_tu[i, j]  <- test[4]
    n_under[i, j]   <- test[5]
    n_over[i, j]    <- test[6]
    
  }
  print(j)
}

ps$place <- 1:183
ps <- melt(ps, id = "place")

p_five <- seq(10,50, 10)

mockvars.labels   <- c(par.nat = "Parents is \nnative",
                       par.fem = "Parent is \nfemale",
                       par.edu = "Parent has \nbachelor's",
                       chi.fem = "Child is \nfemale",
                       chi.old = "Child is \noldest",
                       chi.coh = "Parent and child \ncohabit")  

pcurves <- 
  ggplot(ps, 
         aes(x = place, 
             y = value)) + 
  geom_line() + 
  facet_grid(variable~.,
             labeller = labeller(variable = mockvars.labels)) + 
  geom_hline(yintercept = 0.15, 
             linetype = "dotted") +
  geom_hline(yintercept = 0.05, 
             linetype = "dashed") +
  theme_classic() +
  xlab("Bandwidth in days") + 
  ylab("One-sided p-value")  +
  geom_line(data = data_frame(x=c(1,1,1), y=c(0,0,0), type = c("solid", "dotted", "dashed")),
            aes(x, y, linetype = type)) +
  theme(axis.title = element_text(size = 14),
        legend.position = "top") +
  scale_linetype_discrete(name = "",
                          labels = c("P-values", "p = 0.15", "p = 0.05"))

pcurves_10days <- 
  ggplot(ps %>% filter(place < 11), 
         aes(x = place, 
             y = value)) + 
  geom_line() + 
  facet_grid(variable~.,
             labeller = labeller(variable = mockvars.labels)) + 
  geom_hline(yintercept = 0.15, 
             linetype = "dotted") +
  geom_hline(yintercept = 0.05, 
             linetype = "dashed") +
  theme_classic() +
  xlab("Bandwidth in days") + 
  ylab("One-sided p-value")  +
  geom_line(data = data_frame(x=c(1,1,1), y=c(0,0,0), type = c("solid", "dotted", "dashed")),
            aes(x, y, linetype = type)) +
  theme(axis.title = element_text(size = 14),
        legend.position = "top") +
  scale_linetype_discrete(name = "",
                          labels = c("P-values", "p = 0.15", "p = 0.05")) + 
  scale_x_continuous(labels = 1:10, breaks = 1:10)

